function x = GaussSeidelSolver(lhs, rhs, numIter)    
    x = zeros(size(lhs, 2), 1);
    L = sparse(tril(lhs, 0));
    U = sparse(triu(lhs, 1));
      
    for i = 1 : numIter        
        x = L \ (rhs - U * x);
        residual = lhs * x - rhs;
        mse = residual' * residual;
       %fprintf(1, 'Iteration %d residual %.3f\n', i, mse);        
    end
end